/*
Pittsburgh May 2002
CALCULATE BOND LENGTHS FOR ATOM PAIRS, CONDITIONALLY CONVOLUTE WITH GAUSSIAN
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
typedef double V3[3];
typedef double V33[3][3];

/* XYZ */
#define mxspec 111
#define dim_xyz 13000 /* num_of_at dimension */
#define dim_spec 10 /* num_of_spec dimension */
#define dim_spec_pairs dim_spec*(dim_spec-1)/2
V3 r[dim_xyz];
int ich[dim_xyz],site[dim_xyz];
int spec[dim_spec],nspec[dim_spec];
int specpo[mxspec]; /* pointer to chemistry label */
double nnc[dim_spec_pairs],cnorm[dim_spec_pairs];
int ch2spec[dim_spec][dim_spec];
int ch2ch[dim_spec_pairs][2];
int ext[3];
int nat; /* number of atoms */
int nnspec,nn2spec;
V3 ba[3];

/* GAUSS */
#define nsp 50000
#define nga 1000
double rdf[dim_spec_pairs][nsp];
double gau[2*nga+1];
double rpix;
int npix,ngpix;
int dogauss;

/* cccccccccccccccccccccccccccccccccccccccccccccccccccccccc */

double bound(x)double x;
{double y;
  y=x-floor(x);
  if (y<-0.49999)y=y+1;
  if (y>0.50001)y=y-1;
  return y;
}
double volume(a)V33 a;
{double det;
  det=a[0][0]*a[1][1]*a[2][2]+a[0][1]*a[1][2]*a[2][0]+
      a[1][0]*a[2][1]*a[0][2]-a[0][2]*a[1][1]*a[2][0]-
      a[0][0]*a[1][2]*a[2][1]-a[2][2]*a[0][1]*a[1][0];
  return fabs(det);
}
void pairs(i,j,rc) long i,j; double rc;
{
  long k,o,m,n,p1,p2,p3,icol;
  double bondl,wt;
  V3 w,wr;
  icol=ch2spec[specpo[ich[i]]][specpo[ich[j]]];
  for (p1 = -ext[0]; p1 <= ext[0]; p1++) {
    for (p2 = -ext[1]; p2 <= ext[1]; p2++) {
      for (p3 = -ext[2]; p3 <= ext[2]; p3++) {
	if (i!=j||(i==j&&labs(p1)+labs(p2)+labs(p3)>0)) {
	  for (o = 0; o <= 2; o++){wr[o]=bound(r[j][o]-r[i][o]);}
	  for (o = 0; o <= 2; o++)
	    w[o]=(wr[0]+p1)*ba[o][0]+(wr[1]+p2)*ba[o][1]+(wr[2]+p3)*ba[o][2];
	  bondl=sqrt(w[0]*w[0]+w[1]*w[1]+w[2]*w[2]);
	  if (bondl<rc){
	    if(dogauss==0){
	      printf("%8.4f %8.4f %8.4f %8.4f %4i %4i %3i %3i %3i %3i\n",bondl,w[0],w[1],w[2],i+1,j+1,ich[i],ich[j],site[i],site[j]);
	      if(i!=j){
/*export also this, because later max-gap criterion may break this symmetry*/
		printf("%8.4f %8.4f %8.4f %8.4f %4i %4i %3i %3i %3i %3i\n",bondl,-w[0],-w[1],-w[2],j+1,i+1,ich[j],ich[i],site[j],site[i]);
	      }
	    }else{
	      k=floor(bondl/rpix+0.5);
/*  	      printf("icol ind bondl %i %i %6.2f\n",icol,k,bondl); */
	      wt=2.0; if(i==j){wt=1.0;}
	      for(m=-ngpix;m<=ngpix;m++){
		if((k+m)>=0&&(k+m)<=npix){rdf[icol][k+m]+=wt*gau[m+ngpix];
		}
	      }
	    } /* end print or dogauss */
	  } /* end bondl<rc */
	}
      }
    }
  }
}

main(int argc, char *argv[])
{
  FILE *txt;
  long i,j,m,n,o,p,q;
  char *cln;
  int c,k,ie;
  double cutoff=10.0;
  double basl;
  /* GAUSS */
  int donorm=0;
  int colmax;
  double pi,cg,cc,rr;
  double tmp,sum,vol;
  V3 rtpm;
  double sg=0.1; /* this is sigma */
  double res=0.1; /* bin size in sigma units */

  /* COMMAND LINE ARGS */
  ie=1;
  while(--argc>0){
    cln=*++argv;
    /* NOTE: strcmp is zero on match! */
    if(!strcmp(cln,"-rc")){
      cln=*++argv;
      cutoff=atof(cln);
      --argc;
    }
    if(!strcmp(cln,"-sg")){
      cln=*++argv;
      sg=atof(cln);
      --argc;
      dogauss=1;
    }
    if(!strcmp(cln,"-norm")){ /* normalize */
      donorm=1;
    }
  }

  for(i=0;i<mxspec;i++){specpo[i]=0;} /* init pointers to spec labels */

  /* READ ATOMS: assuming filtered through: xyzfilter -bonds -xyz */
  scanf("%ld%*[^\n]",&nat); /* atoms */
  scanf("%ld%*[^\n]",&nnspec); /* "nnspec" species defined */
  for(i=0;i<nnspec;i++){ /* each "spec" has "nspec" atoms */
    scanf("%ld%ld%ld%*[^\n]",&p,&spec[i],&nspec[i]);
    specpo[spec[i]]=i; /* note this DOES reserve spec=0 for vacancy */
  }
  scanf("%ld%*[^\n]",&nn2spec); /* "nn2spec" species pairs of species */
  for(i=0;i<nn2spec;i++){ /* ch2spec contains pair of species */
    scanf("%ld%ld%ld%ld%*[^\n]",&p,&q,&m,&n);
    ch2spec[specpo[m]][specpo[n]]=i;
    ch2spec[specpo[n]][specpo[m]]=i;
    ch2ch[i][0]=m; ch2ch[i][1]=n;
    if(m!=n){nnc[i]=2*nspec[specpo[m]]*nspec[specpo[n]];}
    else{nnc[i]=nspec[specpo[m]]*nspec[specpo[m]];}
  }
  ext[0]=0; ext[1]=0; ext[2]=0;
  for (i=0; i<3; i++) { /*basis vectors*/
    basl=0;
    for (j=0; j<3; j++){scanf("%lg", &ba[j][i]);basl+=ba[j][i]*ba[j][i];}
    basl=sqrt(basl); if(cutoff>basl/2){ext[i]=floor(cutoff/basl)+1;}
    scanf("%*[^\n]");
  }
  for(i=0;i<nat;i++){ /* read XYZ CHEM SITE */
    for (q=0; q<3; q++){scanf("%lg", &r[i][q]);}
    scanf("%ld%ld%*[^\n]",&p,&q);
    if(p>0){ich[i]=p; site[i]=q;} /*species*/
  }

/*    printf("%i atoms\n",nat); */
/*    for(i=0;i<nnspec;i++){ */
/*      printf("%i %i  spec nspec\n",spec[i],nspec[i]); */
/*    } */
/*    for(i=0;i<nn2spec;i++){ */
/*      printf("%i %i %i  pair\n",i,ch2ch[i][0],ch2ch[i][1]); */
/*    } */
/*    exit(0); */

  /* Conditional gauss init */
  if(dogauss==1){
    pi=4*atan(1.0);
    rpix=sg*res;
    npix=floor(cutoff/rpix+0.5);
    if(npix>nsp){printf("Current resolution res=%1.4f in sigma=%1.4f units and cutoff=%1.1f needs %i pixels, but dimension nsp=%i.\n",res,sg,cutoff,npix,nsp);exit(0);}
    ngpix=floor(5*sg/rpix)+1;
    if(ngpix>nga){printf("Current gauss-dim is %i pixels, need gau=%i.\n",nga,ngpix);exit(0);}
    colmax=nn2spec;
    for(i=0;i<=npix;i++){
      for(j=0;j<dim_spec_pairs;j++){
	rdf[j][i]=0;
      }
    }
    cc=1/sqrt(2*pi)/sg;
    sum=0;
    for(i=0;i<2*ngpix+1;i++){
      j=i-ngpix; cg=j*rpix/sg;
      gau[i]=cc*exp(-cg*cg/2);
      sum+=gau[i]*rpix;
    };
/*      printf("GAUSS: colmax npix sg rpix ngpix %i %i %6.4f %6.4f %i\n",colmax,npix,sg,rpix,ngpix); */
  };

  /* COMPUTE PAIRS: add gauss, or print raw */
  for (i=0; i<=nat-1; i++) {
    for (j=i; j<=nat-1; j++) pairs(i,j,cutoff+2*sg);
  }

  /* conditionally output gauss-convolved data */
  if(dogauss==1){
    vol=volume(ba);
/*      printf("vol rho  %6.1f %7.4f\n",vol,nat/vol); */
    for(j=0;j<colmax;j++){
      cnorm[j]=vol/(4*pi*nnc[j]);
    }
    printf("# %1.2f %1.4f %1.4f    # cutoff sigma resolution\n",cutoff,sg,res);
    printf("# R[A] ");
    for(j=0;j<colmax;j++){printf("%i-%i ",ch2ch[j][0],ch2ch[j][1]);};
    printf("\n");
    for(i=1;i<=npix;i++){
      rr=i*rpix;
      printf("%1.4f ",rr);
      for(j=0;j<colmax;j++){
	tmp=rdf[j][i];
	if(donorm==1){tmp*=1/rr/rr*cnorm[j];}
	printf("%1.6f ",tmp);
      };
      printf("\n");
    }
  }
}
